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1  Introduction 


A  truly  detailed  model  of  an  entire  human  physiological  system  (circulatory,  respiratory,  etc.)  is  currently  not 
feasible  due  to  the  complicated  and  dynamic  geometry,  multi-species  and  multiphase  chemistry,  and  many 
complex  biological  processes  on  a  wide  range  of  spatial  and  temporal  scales.  To  date,  the  main  approach  has 
been  to  create  full  scale,  multi-dimensional,  and  detailed  models  of  a  limited  section  of  a  physiological  system, 
such  as  the  heart,  a  heart  valve,  or  a  section  of  vein  or  artery.  These  models  obtain  the  high-fidelity  solutions 
for  a  small  part  of  the  overall  human  system  and  help  with  local  physiological  diagnoses  and  treatments. 
Another  approach,  taken  here,  is  to  develop  a  broader  and  more  integrated  model  that  facilitates  analysis  of 
larger-scale  material  transport  and  response  in  the  body  and  the  diagnosis  of  time-dependent  systems-level 
phenomena  such  as  shock.  This  approach  also  provides  a  framework  to  integrate  the  calibrated  submodels 
derived  from  detailed  3D  simulations. 

It  is  with  these  goals  in  mind  that  we  develop  a  global,  time  dependent,  yet  low-dimensional  fluid 
dynamics  model  of  the  human  circulatory  system.  Development  of  this  systems-level  model  is  our  first  step 
of  a  larger  effort  to  develop  models  of  individual  physiological  systems  and  then  to  couple  them  together, 
while  also  permitting  the  use  of  both  experimental  data  and  more  detailed  simulations  for  calibration  (Green 
et  al.,  21-23  November,  2010;  Staples  et  ai,  18-20  November,  2007).  Each  individual  model  will  also  allow 
for  future  extensions  that  can  directly  incorporate  more  branching  networks  or  more  detailed  representations 
of  the  chemical  or  physical  processes  occurring.  The  goal  of  this  larger  effort  is  not  to  simulate  details  of 
individual  systems,  but  to  study  the  macro-system  dynamics  of  the  human  body,  as  the  effects  of  coupled 
interactions  may  be  more  important  than  the  details  of  the  individual  physiological  systems.  We  look  to 
determine  the  necessary  level  of  detail  and  accuracy  in  the  submodels  of  an  entire  physiological  system 
needed  to  reproduce  observed  behavior  in  realistic  timeframes. 

Most  physiological  systems  in  the  human  body,  such  as  the  circulatory,  respiratory,  and  lymphatic 
systems,  involve  the  transport  of  gases  and  liquids  through  elastic  volumes  and  channels.  They  are,  therefore, 
governed  by  the  equations  of  reactive  fluid  dynamics.  As  a  first  step,  we  will  consider  a  simplified  model  of 
the  circulatory  system,  shown  schematically  in  figure  1(a),  as  fluid  circuits  that  incorporate  blood  flow  in 
Manuscript  approved  August  22,  2011. 
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systemic  circulation 
(upper  torso,  head,  neck) 


Figure  1:  (a)  Simplification  of  the  circulatory  system,  (b)  Diagram  of  our  initial  representation  of  the  human 
circulatory  system. 


idealized  flexible  channels  and  chambers.  The  system  shown  in  figure  1(b)  is  low-dimensional:  three  spatial 
dimensions  are  collapsed  to  one  dimension  with  elastic  cylindrical  geometry.  The  description  is  dynamic  and 
includes  higher-dimensional  effects  such  as  boundary  layers  and  the  transition  from  convective  to  diffusive 
transport.  This  papers  describes  the  relevant  equations,  the  numerical  schemes  used  to  obtain  fast  and 
efficient  simulations  on  a  laptop  computer,  and  benchmarks  of  the  human  circulatory  system  against  which 
the  model  was  compared. 


2  Background 

Many  numerical  models  attempt  to  recreate  the  human  circulatory  system,  ranging  from  detailed  fluid- 
chemical  simulations  to  lower-dimensional  lumped  parameter  models.  Our  efforts  parallel  those  that  simulate 
biofluid  flows  with  low-dimensional  models.  Among  the  first  was  Womersley,  who  studied  the  propagation 
of  pressure  waves  in  elastic  cylinders  as  a  model  of  arterial  flow  (Womersley,  1955,  1957).  This  work  showed 
that  for  low  forcing  frequency,  small  radius,  low  viscosity,  and  low  elasticity,  it  was  adequate  to  assume  that 
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the  velocity  profile  was  axisymmetric  and  parabolic.  Since  then,  many  efforts  have  used  one-dimensional 
fluid  equations  to  model  the  flow  through  large  sections  of  the  arterial  system  (Schaaf  &  Abbrecht,  1972; 
Zagzoule  &  Marc-Vergnes,  1986;  Sheng  et  al.,  1995).  In  these  cases,  connected  arterial  segments  were  defined 
with  assigned  length,  radius,  and  compliance.  The  one-dimensional  flow  through  each  of  these  segments  was 
calculated  using  boundary  conditions  at  branching  locations  that  were  determined  using  two  principles:  a 
conservation  of  mass  as  the  vessels  branched  and  pressure  equivalence  across  the  branching  point. 

The  complexity  of  these  models  increases  rapidly  with  the  addition  of  more  and  smaller  segments  of 
the  arterial  tree.  Olufsen  et  al.  (2000)  and  Olufsen  (1999)  handled  this  by  first  calculating  the  flow  through 
large  arteries  using  the  nonlinear  Navier-Stokes  equations,  and  then  combining  segments  of  smaller  arteries 
together  as  a  structured  tree  model  governed  by  the  linearized  equations.  This  approach  generated  appropri¬ 
ate  outflow  boundary  conditions  for  the  nonlinear  calculations  in  the  large  arteries,  and  the  results  compared 
favorably  with  experimental  measurements. 

Similarly,  low-dimensional  models  of  larger  arterial  geometries  provided  boundary  conditions  for  more 
highly  resolved  computational  simulations  (Vignon-Clementel  et  al.,  2006).  This  approach  also  showed 
flow  rates  computed  by  the  one-dimensional  vascular  models  that  compared  well  with  fully  resolved  three- 
dimensional  simulations  and  in  vivo  measurements  (Steele  et  al.,  2003;  Wan  et  al.,  2002).  From  the  success 
of  these  results,  the  authors  speculated  that  such  low-dimensional  models  could  be  of  direct  medical  use 
by  providing  a  tool  for  real  time  surgery  and  therapy  planning.  Grinberg  et  al.  (2009)  performed  a  three- 
dimensional  simulation  of  the  circulation  in  the  intracranial  arterial  system,  and  were  also  able  to  show  good 
agreement  in  distributions  of  pressure  and  mass  flow  when  compared  against  a  simpler  one-dimensional 
model  (Grinberg  et  al.,  2011). 

Low  dimensional  models  have  also  been  constructed  using  electrical  circuit  components,  such  as  capac¬ 
itors  and  resistors,  to  simulate  the  circulatory  system  dynamics  (Noordergraaf  et  al.,  1963;  Snyder  et  al., 
1968;  Westerhof  et  al.,  1969).  These  circuit  models  used  current  to  represent  one-dimensional  fluid  flow  and 
electric  potential  to  represent  the  pressure  differential  along  vessels.  Electrical  circuits  were  physically  con¬ 
structed  as  an  “analog  computer”  that  produced  voltage  and  current  data  that  compared  well  with  measured 
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blood  pressure  and  flow  rates. 


As  an  example  of  a  comprehensive  physiological  that  is  not  essentially  based  on  the  equations  of  fluid 
dynamics  is  the  Physiome  Project  at  the  National  Simulation  Resource  (NSR).  The  Physiome  project  is 
a  collection  chemical  and  biological  models  at  a  range  of  levels  and  an  international  effort  to  use  them  to 
create  a  full  description  of  the  human  physiome  (Neil  &  Bassingthwaighte,  2007;  Bassingthwaighte,  2000). 
NSR  has  focused  on  creating  open-source,  user-friendly  software  and  standard  data  formats,  all  of  which 
encourages  collaboration  and  modularization  of  the  models  available  there. 

Our  circulatory  system  model  is  a  closed  circuit  in  which  axisymmetric  flow  is  pumped  by  a  periodic 
contraction  and  relaxation  of  the  heart  chamber  geometry.  Prior  models  drove  the  fluid  flow  using  prescribed 
boundary  conditions  on  the  pressure  or  velocity.  Sheng  et  al.  (1995)  for  example,  provided  the  pulsatile 
aortic  and  constant  caval  (venous)  pressures  as  inlet  and  outlet  conditions  of  their  model  system.  Other 
physiological  phenomena,  such  as  the  implementation  of  valves  to  inhibit  backflow  in  the  heart  chambers  and 
the  low-pressure  venous  system,  and  augmented  blood  pumping  provided  by  the  skeletal  muscular  system 
throughout  the  body,  are  also  included  in  the  development  of  our  model  circulatory  system  and  differentiates 
our  work  from  previous  efforts. 


3  The  physical  model 

A  diagram  of  the  model  circulatory  system,  figure  1(b),  shows  the  four  chambers  of  the  heart:  the  right 
atrium  and  right  ventricle  (RA,  RV),  and  the  left  atrium  and  left  ventricle  (LA,  LV).  There  are  four  heart 
valves  (mitral,  aortic,  tricuspid,  and  pulmonary)  that  control  the  flow  through  the  heart  and  into  the  systemic 
and  pulmonary  systems.  Semi-lunar  valves  in  the  vasculature  prevent  flow  reversal  which  might  occur  due  to 
the  pulsatile  nature  of  the  heart  pumping  and  resulting  drops  in  pressure.  The  pulmonary  and  vascular  trees, 
shown  schematically  on  the  top  and  bottom  of  the  figure,  are  comprised  of  elastic  vessels  and  contain  much  of 
the  interesting  physics  of  the  problem.  Here,  the  system  branches  form  the  large  arteries  into  smaller  arteries, 
into  the  smaller  arterioles,  and  finally  to  the  capillaries.  To  return  to  the  heart,  the  capillary  beds  coalesce 
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into  the  ventricles,  collect  into  the  smaller  veins,  and  deposit  into  the  large  veins  (inferior  and  superior  vena 
cava).  As  the  fluid  moves  back  and  forth  between  large  vessels  and  smaller  capillaries,  the  inherent  flow 
physics  change  from  being  dominated  by  convective  effects  to  being  dominated  by  viscous  effects. 

Beginning  at  the  left  atrium  of  the  heart  (LA),  oxygenated  blood  (red)  is  driven  through  the  mitral 
valve  into  the  left  ventricle  (LV),  and  then  through  the  aortic  valve  into  the  systemic  section  of  the  circuit. 
When  the  flow  enters  the  systemic  vascular  tree,  the  total  cross  sectional-area  increases  and  there  is  a  large 
pressure  drop  to  drive  the  flow  through  the  capillaries,  where  oxygen  is  delivered  to  organs,  muscles,  etc. 
After  passing  through  the  systemic  vascular  tree,  the  deoxygenated  blood  (blue)  then  passes  through  the 
venous  system,  and  is  drawn  into  into  the  right  atrium.  From  here  it  is  pumped  through  the  tricuspid 
valve  into  the  right  ventricle,  then  through  the  pulmonary  valve  into  the  pulmonary  vasculature,  where  a 
smaller  yet  comparable  pressure  drop  drives  the  flow.  Here,  the  pulmonary  capillaries  are  intertwined  with 
the  small-scale  alveoli  of  the  pulmonary  system,  exchanging  oxygen  and  carbon  dioxide  with  the  lungs.  The 
oxygenated  blood  (red)  is  returned  to  the  left  atrium  via  the  pulmonary  veins,  completing  the  circuit. 


3.1  Flexible  Flow  Equations 


The  mathematical  formulation  for  the  flow  is  based  on  a  solution  of  the  Flexible  Flow  Equations  (FFEs). 
Given  appropriate  initial  and  input  conditions,  the  equations  are  solved  for  the  pressure,  velocity,  and  cross- 
sectional  area  of  an  incompressible  flow  as  a  function  of  time  and  position  along  a  system  of  elastic  channels. 
By  assuming  an  axisymmetric  parabolic  velocity  profile,  the  physical  system  collapses  from  three  spatial 
dimensions  and  time  to  one  spatial  dimension  and  time. 

The  fluid  velocity  u{S,  t)  at  a  distance  S  along  the  loop  and  time  t  is  described  by  a  momentum  equation. 


du{S,  t) 
dt 


u{S,  t) 


du{S,  t) 
dS 


1  dp{S,  t)  d'^u{S,  t) 

p  dS  dS^ 


(1) 


where  p  is  the  density,  v  is  kinematic  viscosity,  R  is  the  local  radius  of  the  pipe,  g  is  an  external  acceleration 
from  gravity  or  any  applied  body  force,  and  p  is  the  local  pressure.  Equation  f  satisfies  the  parabolic  viscous 
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similarity  solution  with  zero  velocity  at  walls. 

Because  the  fluid  velocities  are  low  compared  to  the  speed  of  sound,  we  can  assume  that  the  fluid 
acceleration  is  small  compared  to  the  pressure  gradient.  Equation  1  is  rearranged  as, 


Sniy  1  dp  du  du  d"^u  1  dp  , 


(2) 


where  the  generalized-acceleration  term  G{S,  t)  is  defined  as, 


G{S,t)  =  9 


du 

m 


+  V 


d^u 


(3) 


In  one  limit,  we  can  assume  that  G  =  0,  so  that  this  acceleration  is  truly  small  compared  to  the  pressure- 
gradient  term.  This,  however,  is  generally  not  valid  in  the  cases  we  are  considering,  where  the  acceleration 
can  be  large  during  fluid  pumping.  In  §4  of  this  paper,  when  we  present  a  solution  method  for  the  equations, 
we  describe  a  method  for  including  portions  of  the  acceleration  term.  We  also  present  and  compare  results 
for  cases  with  G  =  0  and  G  7^  0. 

Defining  A  =  and  rearranging  equation  2  gives, 

(4) 

^TTV  ^TlVp  dS 


Flexible  walls  are  modeled  by  allowing  the  local  area  ^(S',  t)  to  accommodate  the  local  pressure  through 
an  equation  of  state, 

P~P<^<i  =  ~  •  (5) 

Here,  Agq  is  the  cross-sectional  area  of  the  the  flexible  channel  at  the  equilibrium  pressure  Peq,  and  s{S)  is 
a  dimensionless  elasticity  coefficient  defined  along  the  length  of  the  channel.  As  the  elasticity  goes  to  zero 
(rigid  walls),  the  area  is  pegged  at  the  equilibrium  area  Aeq,  independent  of  pressure  changes.  Note  that  Aeq 
may  be  a  function  of  S  and  t,  and  is  used  in  the  model  of  the  vascular  trees  as  discussed  below. 
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Even  though  the  flow  is  incompressible,  the  local  mass  flow  rate  is  not  necessarily  constant  because  the 
walls  are  elastic.  The  total  mass  of  the  fluid  in  the  system  remains  constant,  however,  and  therefore  any 
local  expansion  must  be  balanced  by  contraction  elsewhere.  This  is  reflected  in  the  continuity  equation. 


dA 


(6) 


Combining  equations  5  and  6  yields  the  pressure  equation, 

^  _  Peg  dA  _  -Peg  d  {uA) 
dt  sAeg  dt  sAeg  dS 

The  term  uA  can  be  obtained  through  equation  4, 

,  A^G  A^  dp 

qiA  —  _  _ _ — 

SiTiy  Swiyp  dS ' 

Combining  these  equations,  we  arrive  at  a  second-order  diffusion  equation  for  the  pressure, 


dp{S,t)  Peg  d  J 

{  AHS,t)  dp{SA)] 

[  Peg  d  1 

{ AHSA)GiS,t)\ 

dt  ~  e{S)AegiS)dS] 

18tti/{S)p  dS  J 

\  e{S)A,giS)  dS  ] 

1  8TTiy{S)  j 

(7) 


(8) 


(9) 


3.2  Convective  to  diffusive  flow  in  vascular  trees 

As  the  blood  vessels  branch  in  the  systemic  and  pulmonary  vascular  trees,  the  total  cross-sectional  area 
increases  rapidly  and  the  flow  rate  drops  correspondingly.  Each  vessel  through  which  blood  flows  decreases 
in  size,  and  therefore  viscous  effects  become  increasingly  important.  To  address  the  transition  from  convective 
to  diffusive  flow  in  the  systemic  and  pulmonary  vascular  trees,  we  define  a  weighted  viscosity  term,  1^(8), 
that  artificially  increases  in  the  small-scale  capillary  regions  where  the  viscous  term  dominates  despite  the 
increasing  A.  The  weighted  viscosity  is  a  smooth  function  that  varies  from  a  value  of  blood  viscosity  at 
the  upstream  and  downstream  ends  of  the  vascular  tree  to  a  maximum  value  at  the  widest  section  of  the 
vasculature  representing  the  capillary  bed,  with  the  maximum  value  calibrated  to  reproduce  a  reasonable 
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Figure  2:  Area  and  viscosity  along  the  length  of  the  circulatory  system. 


pressure  drop  across  the  vascular  tree.  The  weighted  viscosity  is  written  as, 


iy{S)  =  Vb 


[S  -  S{XMS{Xn)  -  Sf\ 
{c{S{Xm)-S{X,)YY  )  ’ 


(10) 


where  is  the  blood  viscosity,  c  is  a  constant  normalizing  the  distribution  term  to  v^ax  at  the  peak,  and 
S'(Ai)  and  S{Xm)  are  the  locations  of  the  upstream  and  downstream  ends  of  a  vascular  tree. 

The  equilibrium  area  is  initialized  using  a  similar  function  to  the  viscosity  in  the  vascular  trees.  The 
variation  of  both  the  equilibrium  area  and  the  viscosity  along  the  length  of  the  circulatory  model  are  shown 
in  figure  2.  In  this  case,  Xi  =  0.5m  and  X^  =  2m  in  for  the  systemic  vascular  tree,  while  Xi  =  3m  and 
Xjs!  =  4.5m  for  the  pulmonary  vascular  tree.  According  to  Guyton  &  Hall  (2000),  about  5  L  is  considered 
normal  for  an  adult  human,  and  the  area  distribution  of  our  current  model  yields  a  total  blood  volume  of  4.7 
L.  The  distribution  of  weighted  viscosity  and  equilibrium  area  across  the  model  circulatory  system  is  shown 
in  figure  2. 


3.3  Heart  and  skeletal  muscle  pumping 

The  blood  flow  is  pumped  by  a  prescribed  periodic  contraction  and  relaxation  of  the  heart  chambers,  and 
the  local  time-varying  pressure  is  not  prescribed  anywhere.  The  effects  of  muscle  contraction  and  relaxation 


on  the  circulatory  system  are  represented  through  imposed  local  changes  in  the  equilibrium  area.  At  each 


ventricle 


Figure  3:  (a)  Cross-sectional  area  of  the  atria  and  ventricles  during  the  period  of  one  hearbeat.  (b)  Cross- 
sectional  area  of  the  vascular  tree  (0.5  <  S  <  2)  area  during  skeletal  muscle  contraction.  Contraction  only 
applied  on  the  venous  end  of  the  vascular  tree  (1.25  <  5”  <  2).  For  the  purposes  of  the  visualization,  the 
maximum  amplitude  of  contraction  is  AAm  =  0.35,  and  the  frequency  is  1  Hz. 


location  and  time,  Af,q{S,t)  is  either  increased  or  decreased  to  model  a  local  relaxation  or  contraction  of 
a  muscle,  respectively,  and  the  local  pressure  is  then  adjusted  accordingly.  In  the  atria  of  the  heart,  for 
example,  the  contraction  occurs  for  time  {tai  <  t  <  ta2),  where  (ta2  ~  iai  =  2a)  is  the  period  of  atria 
contraction.  Similarly,  Ty  is  the  period  of  ventricle  contraction.  The  corresponding  equilibrium  area  change 
in  the  atria  is  modeled  as. 


Aeq{Sa,t)  =  Aqo{Sa){^  “  cr(t)(l  -  AH^)),  where 


cr(t) 


6.75 


/  {t  -  tai){ta2  - 

V  {tal  -  ta2Y 


1.5 


(11) 


Here,  AAa  is  the  fraction  of  the  maximum  area  change  during  the  atrial  contraction  (AH„  in  the  ventricles). 
The  areas  of  the  atria  and  ventricles  in  the  heart,  over  the  duration  of  one  heartbeat,  are  shown  in  figure  3(a). 


General  systemic  contractions  of  skeletal  muscles  are  represented  similarly.  Small  amplitude  periodic, 


distributed  contractions  are  applied  to  the  venous  side  the  systemic  vascular  tree  (5'(Amid)  <  S  <  S{X]y), 
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where  S{Xmid)  =  0.5[S'(Xi)  +  S'(Xjv)]).  Both  the  magnitude  and  the  frequency  are  adjustable,  so  that 

27rt 

^mt  = 

^  m 

2Tr{S  -  SiXi)){S{XN)  -  S) 

(SiX^)  -  Six,)) 

^eq  —  ^eqO  (l  “t“  COsiiOrnt)  Sin(ca77ia,)^  ,  (1^) 

where  is  the  period  of  muscle  contraction,  occurring  during  time  {tam  <  t  <  tbm),  and  AAm  is  the 
maximum  contraction  amplitude.  A  visualization  of  the  muscle  contraction  model  is  shown  in  figure  3(b). 

Valves  are  an  important  part  of  the  physical  circulatory  system  to  prevent  flow  reversal.  This  is  necessary 
within  the  heart  chambers,  driving  fluid  out  through  the  arteries  instead  of  back  through  the  atria  and  veins. 
In  the  low-pressure  venous  region,  semi-lunar  valves  are  also  integral  in  directing  blood  back  towards  the 
heart  by  preventing  backflow  and  pooling  in  the  extremities.  In  the  current  work,  both  heart  valves  and 
venous  semi-lunar  valves  are  actuated  in  the  same  manner.  At  each  valve  location,  the  local  pressure  is 
monitored  at  every  timestep.  When  an  adverse  pressure  gradient  is  detected,  the  local  area  A  is  initially 
decreased  by  80%,  and  then  to  0%  at  the  next  time  step  if  the  adverse  pressure  gradient  persists. 

4  Numerical  Implementation 

Figure  4  is  a  schematic  of  the  layout  of  the  computational  domain  based  on  Figure  1(b).  Various  regions 
of  the  circulatory  system  can  be  represented  by  an  arbitrarily  chosen  number  of  computational  cells,  thus 
providing  more  or  less  resolution  in  those  regions,  as  discussed  in  §5.  This  initial  representation  is  a  single 
loop,  but  the  algorithm  can  be  extended  to  treat  different  branches  and  sub-loops  in  a  more  complex  network. 

Equations  5,  8,  and  9  are  discretized  to  produce  a  tridiagonal  system  that  can  be  solved  iteratively  at 
each  time  step.  Equation  9  is  discretized  using  a  time  derivative  of  pressure  and  a  central  spatial  derivative 
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MV  AV  SLV  SLV  TV  PV  SLV  SLV 


Figure  4:  Diagram  of  the  numerical  implementation  of  the  model.  At  top:  LA,  left  atrium;  MV,  mitral  valve; 
LV,  left  ventricle;  AV,  aortic  valve;  SLV,  semi-lunar  valves;  RA,  right  atrium;  TV,  tricuspid  valve;  RV,  right 
ventricle,  and  PV,  pulmonary  valve.  Top  scale  line  gives  the  dimensional  distances,  and  the  bottom  scale 
line  shows  the  distribution  of  computational  cells  for  the  baseline  case,  described  in  §5.  Three  cell-centered 
green  targets  indicate  the  approximate  locations  in  the  arteries  (SA),  the  capillaries  (SC),  and  the  veins 
(SV),  at  which  data  is  sampled  in  §5. 


of  a  general  variable  /, 


dp{S,  t) 
dt 


1 

5t 


5/ 

dSj 


5S. 


(13) 


where  /  can  be  any  of  the  model  quantities  that  vary  in  space  (/(S')).  The  subscript  j  indicates  spatial 
indices  in  the  range  [0,  /iV-|-l],  where  N  is  the  number  of  computational  cells,  and  the  superscripts  k  and 
k—1  indicate  the  current  and  previous  timestep,  respectively.  The  areas  and  velocities  are  calculated  at  cell 
interfaces  and  at  the  half  timestep,  so  that  the  time  rate  of  change  of  pressure  can  be  solved  as  a  central 
difference. 

Equation  9  then  becomes 


+Pj{^  +  +  i  -  ^J^J  +  hPj  +  1  =  Pj 


k-1 


(14) 


where 


8TTiy^+ip6Sj+i  ’ 


a  ^  St  Peg 


^3  +  k  = 


8-nv4 


'l+3 


(15) 
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Equation  14  is  arranged  into  a  tridiagonal  system, 


(16) 


with  vector  coefficients, 


Aj  =  —/3jaj_i 

Bj  =  1  +  +  a^-i ) 


(17) 


This  matrix  system  is  solved  using  a  fast  tridiagonal  solver  (Boris,  1976),  and  areas  and  velocities  are 
calculated  using. 


A^  = 


{Pj  -  Peg) 


Peg 


Ak 


STTl/jpSSj 


(18) 


The  procedure  is  iterated  until  the  pressure  solution  converges. 


Including  the  effects  of  the  generalized  acceleration  term  G  complicates  the  solution  process  significantly. 
Assume  for  now  that  the  spatial  gradients  are  small  compared  to  the  temporal  gradients,  and  that  there  are 


no  external  body  forces.  Then 


(19) 


Applying  this  to  equation  4,  yields. 


A  du  A  dp 

Stti'  dt  Sttvp  dS 


(20) 
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Discretizing  equation  20  yields, 


^+1  =- 


S'!rVj_^_i  St 


SiTiy^+ipSS^+i 


(21) 


and  therefore,  u*  can  be  expressed  as. 


D+ 


1  + 


A>‘  , 

STTl/St 


^3  +  1^  3  +  h  J+ _ 

87r:z^-+i(5t  87ri^^-+ip(5S'^-+i 


(22) 


We  discretize  equation  7  to  obtain. 


Pi-p',-'  = 


-PeqSt 

^-^eqSSj 


(«7l)^;|  -  (nA)J:| 


(23) 


Substituting  equation  22  into  23  yields. 


P3 


/-i 


^3+h-^-^3  +  h  iP3  +  l-P3) 


+  P3VJ-I 


7,-1  ■ 


_£ _ 2_ 

(5t 


iP3-P3-l) 


,  (24) 


where. 


«i+i  = 


5  87r:y^-+ip(5^j-+i 


/3.,= 


Stpe 


87rK 


3+2 


^3  +  . 


=  1  + 


^±1 

SiTiySt 


-1 


(25) 


Rearranging  terms  yields 


-pjTlj.ia.i 


p1-i  + 


1  +  Pj 


^3  +  h°^3+h-^3-h^3-h 


(26) 


13 


Now,  the  set  of  vector  coefficients  have  become, 


■^3  ~ 

B,  =  1  +/3j(r7j+iaj+i  +V^_ia^_i) 


^3  — 


D,  =  - 


■'3  1  %+h^3  +  k^^ 


Vi-ili-h 


(27) 


which  are  used  in  the  tridiagonal  system  in  equation  16. 

This  generalization  of  the  method  to  G  yf  0  complicates  the  vector  coefficients,  but  provides  the  ability 
to  include  body  forces  and  temporal  terms.  Including  large  spatial  gradients,  however,  requires  further 
development  and  testing,  and  so  is  left  for  future  research.  A  summary  of  the  numerical  method  is  given  in 
table  1. 


5  Results 

In  this  section,  we  Hrst  present  results  for  a  baseline  simulation  that  we  use  for  comparison  throughout  the 
section.  This  simulation  uses  the  simplest  physical  model,  which  assumes  that  the  acceleration  term  is  equal 
to  zero  (G  =  0).  The  simulation  was  also  performed  using  the  G  =  —du/dt  acceleration  model,  and  they  are 
compared  to  the  G  =  0  case  in  §5.2.  Additionally,  we  present  the  model  response  to  varying  vessel  elasticity 
(§5.3)  and  application  of  the  generalized  muscle  contractions  (§5.4). 

For  the  cases  presented  here,  the  computational  run  times  were  on  the  order  of  100  times  faster  than 
realtime  (neglecting  time  needed  for  data  output).  Using  the  zero-acceleration  model  and  a  heartbeat  of  Is, 
12  heartbeats  were  computed  in  a  CPU  time  of  60ms  (5ms  per  heartbeat).  When  the  nonzero-acceleration 
model  was  employed,  140  ms  of  CPU  time  were  needed  to  compute  12  heartbeats  (12  ms  per  heartbeat). 
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1.  Initialize  variables  Vj,ej,Aeqj,  Aj,  Uj  at  cell  centers  j. 

2.  Start  time  step  loop.  Current  variables  denoted  by  superscript  k. 

(a)  Define  variable  at  previous  timestep  k  —  1,  =  p^ =  A^) 

(b)  Apply  equations  11  and  12  at  locations  of  muscle  contractions 

•  Prescribed  A'lq  j  in  contracting  cells,  resultant  calculated  using  equation  5 

(c)  Activate  valves  if  needed 

•  At  prescribed  locations,  check  for  adverse  pressure  gradient  >  {pjv-i  +  10.0)) 

•  If  so,  decrease  local  area  by  80% 

•  If  already  decreased  in  previous  timestep,  set  local  area  to  0 

(d)  Activate  distributed  muscle  contractions 

•  Prescribed  A^^  j  in  user-defined  region,  resulting  pj  calculated  using  equation  5 

(e)  Solve  for  p  distribution 

•  Calculate  variables  at  interfaces  {A^._  i ,  u^._  i ) 

^2  3  2 

•  Set  up  tridiagonal  system  using  equations  in  §4 

•  Calculate  p^  using  tridiagonal  solver.  Update  A^  and  Uj  using  equation  18 

(f)  If  \pj  —  >  ptoi  (user-defined  tolerance)  anywhere,  update  =  p^  and  returning 

to  step  5. 

(g)  If  \pj  —  pf^^\  <  ptoi  everywhere,  end  time  loop.  Return  to  2. 


Table  1:  Numerical  method  summary 
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5.1  G  =  0 


The  baseline  simulation  is  initialized  at  t  =  0  with  zero  velocity  and  an  equilibrium  pressure  of  1.01  Patm-, 
which  corresponds  to  a  gauge  pressure  (above  atmospheric)  of  7.6  mmHg.  Blood  kinematic  viscosity  and 
density  are  taken  as  0.05  cm^/s  and  1.06  g/cm^,  respectively.  Figure  2  shows  the  equilibrium  area  and 
weighted  viscosity  along  the  length  of  the  system.  The  elasticity  e  is  determined  empirically  and  set  to  1.0. 
The  distributed  skeletal  muscle  contraction  model  is  not  employed. 

This  simulation  uses  100  computational  cells  to  resolve  a  circulatory  system  that  is  5  m  in  length  with  a 
total  volume  of  4.7  L,  considered  normal  for  an  average  adult  human.  The  distribution  of  these  cells  is  shown 
in  figure  4.  Each  vascular  tree  uses  30  cells  and  begin  six  cells  downstream  of  their  respective  ventricles  (left 
ventricle  for  the  systemic  tree,  and  right  ventricle  for  the  pulmonary  tree).  The  atria  are  located  thirteen 
cells  downstream  of  the  vascular  trees.  Each  chamber  of  the  heart  uses  only  one  computational  cell.  Each 
of  the  four  heart  valves  (mitral,  aortic,  tricuspid,  and  pulmonary)  exists  at  the  downstream  interfaces  of  the 
atria  and  ventricles.  There  are  four  venous  semi-lunar  valves,  two  in  each  half  of  the  circulatory  loop.  One 
valve  is  located  two  cells  downstream  of  the  vascular  tree,  and  one  is  located  2  cells  upstream  of  the  atria. 

Three  sensors  have  been  “inserted”  at  selected  locations  to  describe  the  behavior  of  the  arteries,  capil¬ 
laries,  and  veins.  The  approximate  positions  of  the  sensors  are  shown  as  green  targets  in  figure  4.  Arterial 
data  is  obtained  at  S'  =  0.3  m,  which  is  two  computational  cells  downstream  of  the  right  ventricle,  and  four 
cells  upstream  of  the  systemic  vascular  tree.  The  capillary  data  sensor  is  positioned  at  the  point  of  largest 
area  expansion  (S(Xmid)  =  1-25  m),  and  the  venous  data  sensor  is  located  5  computational  cells  downstream 
of  the  systemic  vascular  tree  at  S  =  2.25  m. 

Profiles  of  the  varying  cross-sectional  area,  gauge  pressure,  and  velocity  at  the  three  sensor  locations 
are  shown  in  figures  5(a)  and  5(b).  They  are  shown  for  the  very  first  cycle  of  the  simulation.  When  the 
data  were  phase-averaged  over  the  first  25  heartbeats,  the  largest  standard  deviation  in  pressure,  among  all 
locations,  was  only  0.0057.  Area  and  pressure  profiles  are  qualitatively  very  similar  to  each  other  because 
these  quantities  are  linearly  related  through  the  equation  of  state  (equation  5).  The  pressure  at  all  three 
locations  begins  at  equilibrium  (7.6  mmHg). 
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(a)  (b) 

Figure  5:  Time  resolved  plots  of  (a)  pressure,  area,  and  (b)  velocity  at  three  locations  in  the  vasculature. 
G  =  0. 


Figure  3(a)  shows  that  the  atria  begin  to  contract  at  0.15s,  causing  the  arterial  pressure  and  area  to 
increase  quickly  (solid  line).  The  pressure  in  the  capillaries  (dashed)  and  veins  (dash-dot)  begin  increasing 
more  slowly.  At  0.25s,  the  ventricles  begin  to  contract,  while  at  0.36s  the  atria  have  begun  to  relax  and 
expand  slowly.  This  atria  expansion  causes  a  brief  drop  in  the  venous  pressure  at  0.36,  but  the  increasing 
ventricle  contraction  increases  the  pressure  until  0.45s.  The  ventricles  begin  to  expand  at  0.44s,  but  the 
closing  of  the  aortic  valve  keeps  the  pressure  from  continuing  to  drop  so  quickly  in  the  arteries.  After  this 
time,  the  expansion  of  both  chambers  causes  a  low  pressure  in  the  heart,  but  the  favorable  gradient  keeps  the 
mitral  and  venous  semi-lunar  valves  open,  and  hence  the  pressure  in  the  veins  drops  to  below  the  equilibrium 
pressure.  It  then  increases  slowly  back  to  Peq  at  0.9s,  after  both  chambers  of  the  heart  have  relaxed  back 
to  the  default  equilibrium  area.  At  the  same  time,  the  pressure  in  the  arteries  and  capillaries  decreases 
smoothly  back  to  the  equilibrium  pressure. 

There  is  a  precipitous  drop  in  peak  pressure  between  the  arteries  and  the  capillaries.  In  the  veins, 
pressure  decreases  further,  and  even  falls  below  the  equilibrium  pressure.  Guyton  &  Hall  (2000)  reports  that 
normal  mean  pressures  are  100  mmHg  in  the  arteries,  20  mmHg  in  the  capillaries,  and  only  5  mmHg  in  the 
veins.  The  curves  presented  here  thus  under-predict  the  arterial  and  capillary  pressures  and  over-predict  the 
venous  pressures.  This  indicates  that  a  more  refined  distribution  of  equilibrium  pressure  and  actuation  of 
the  valve  models  is  needed  and  gives  some  indication  of  future  work. 
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Time-resolved  velocities  are  shown  in  figure  5(b).  Arterial  velocities  (solid  line)  increase  as  the  atria 
begin  to  contract  at  0.15s.  As  the  rate  of  atria  contraction  slows,  the  velocity  in  the  arteries  falls  slightly, 
but  remains  positive.  Velocity  in  the  capillaries  (dashed),  which  is  scaled  using  the  ordinate  on  the  right, 
increases  similarly  at  0.15s,  although  at  a  slower  rate.  Both  velocities  increase  greatly  as  the  ventricle  starts 
to  contract  at  0.25s  Arterial  velocity  peaks  at  0.31s  and  capillary  pressure  peaks  shortly  after  at  0.35s.  As  the 
atria  begin  to  contract,  there  is  a  brief  negative  velocity  in  the  veins  (dash-dot),  but  the  associated  adverse 
pressure  gradient  is  detected  and  mitigated  quickly  by  the  semilunar  valves  in  the  veins.  The  velocity  in  the 
veins  is  relatively  low  during  the  ventricle  contraction. 

As  the  atria  expand  after  0.34s,  the  velocity  in  the  veins  is  small  and  positive,  but  the  venous  velocity 
increases  greatly  at  0.44s,  after  which  all  chambers  of  the  heart  are  expanding.  This  creates  the  low  pressure 
observed  in  figure  5(a),  drawing  blood  through  the  veins  into  the  heart.  In  the  capillaries,  the  velocity 
has  two  peaks,  each  lining  up  with  the  peak  in  the  arterial  velocity,  associated  with  the  contraction  of  the 
heart,  and  the  venous  velocity,  which  occurs  later  with  the  heart  expansion.  In  other  words,  flow  through 
the  capillaries  increases  when  it  is  both  pushed  from  the  arteries  and  pulled  from  the  veins.  Peak  velocities 
in  the  capillaries,  however,  are  much  lower  that  in  the  arteries  and  veins.  This  is  a  result  of  the  increased 
cross-sectional  area  in  the  vascular  trees. 

Time-averaged  velocities  are  compared  with  data  from  Charm  &  Kurland  (1974)  in  table  2  and  average 
values  are  on  the  same  order  of  magnitude.  The  main  discrepancy  is  the  mean  velocities  of  the  arteries 
and  veins.  In  the  reference  data,  the  mean  venous  velocity  is  approximately  half  the  mean  arterial  value, 
while  for  our  baseline  simulation,  the  mean  venous  velocity  is  the  same  or  slightly  larger  than  mean  arterial 
velocity.  Figure  5(b)  shows,  however,  that  the  maximum  velocity  in  the  arteries  is  considerably  larger  than 
the  maximum  venous  velocity.  The  ratio  of  maximum  velocities  (230/129  =  1.8)  is  about  the  same  as  the 
ratio  of  mean  velocities  from  the  reference  data  (45/24  =  1.9). 
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G  =  0 

^  ^max 

G  =  -5u/5t 

U  Umax 

C  K 

u 

Arteries 

30.6  229.7 

49.5  260.8 

45 

Capillaries 

1.2  4.1 

2.0  8.3 

1.7 

Veins 

32.1  128.7 

53.0  173.4 

24 

Table  2:  Mean  velocity  comparison  with  Charm  &  Kurland  (1974).  All  velocities  are  in  cm/s. 
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Figure  6:  Time  resolved  plots  of  (a)  pressure,  area,  and  (b)  velocity  at  three  locations  in  the  vasculature. 
Results  are  from  a  simulation  using  the  acceleration  model  (G  =  —du/dt). 


5.2  Nonzero-acceleration  model 


Using  the  same  input  parameters  and  computational  conhguration,  a  simulation  was  also  calculated  using 
the  nonzero-acceleration  model  described  in  §4.  Here,  the  acceleration  term  G  is  modeled  as  G  =  —du/dt. 
Profiles  of  area,  pressure,  and  velocity  are  obtained  at  the  same  three  “sensor”  locations  and  shown  in 
figure  6.  These  profiles  are  from  the  second  heartbeat  calculated  (Is  <  t  <  2  s).  Phase-averaging  over  24 
heartbeats,  from  the  second  to  the  twenty-fifth,  yielded  a  maxiumum  standard  deviation  in  pressure,  among 
all  locations,  of  only  0.007.  The  first  heartbeat  deviates  from  the  subsequent  cycles  at  the  initial  condition 

(p(0)  =Peq)- 

Peak  pressures  in  the  arteries  are  much  greater  than  those  in  the  capillaries,  which  is  similar  to  the  result 
for  G  =  0.  Also,  the  pressure  in  the  veins  is  less  than  the  equilibrium  pressure  for  most  of  the  heartbeat. 
There  are  significant  differences,  however,  in  both  the  magnitude  of  the  peak  pressures,  and  the  qualitative 
characteristics  of  the  time-varying  profiles. 
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At  the  beginning  of  a  heartbeat,  the  pressures  in  both  the  arteries  (solid  line)  and  the  veins  (dash-dot) 
are  increasing  slightly,  and  a  small  magnitude,  high  frequency  oscillation  is  present  in  the  arterial  pressure. 
Both  of  these  phenomena  are  due  to  the  effects  of  the  previous  heartbeat.  When  the  atria  begin  to  contract 
at  1.15s,  the  arterial  pressure  increases  to  a  second  peak.  As  the  ventricles  begin  to  contract  at  1.25s, 
the  arterial  pressure  increases  more  quickly,  and  reaches  its  maximum  at  1.33s,  the  same  phase  at  which 
peak  arterial  pressure  is  reached  when  using  the  zero-acceleration  model.  The  arterial  pressure  from  the 
acceleration  model  decreases  more  quickly,  however,  dropping  below  the  equilibrium  pressure  at  1.38s  while 
the  ventricle  is  still  contracting.  At  1.34s,  the  atria  begin  to  expand  and  the  pressure  in  the  veins  drops 
below  the  equilibrium  pressure.  The  pressure  in  the  capillaries  (dashed)  increases  slowly  to  its  maximum  at 
1.43s,  at  which  time  the  ventricles  stop  contracting  and  begin  to  expand. 

As  the  ventricles  approach  maximum  contraction,  the  mitral  and  aortic  valves  react  to  the  changing 
pressure  distribution  in  the  heart,  and  their  opening  and  closing  cause  the  fluctuation  in  the  arterial  pressures 
observed  from  1.4s  to  1.43s.  A  longer  time-scale  oscillation  is  observed  in  both  the  arterial  and  venous 
pressures  from  1.45s  until  the  end  of  the  heartbeat  at  2s,  which  is  caused  by  the  inclusion  of  the  acceleration 
model,  as  expressed  in  equation  9.  The  rate  of  change  of  pressure  depends  not  only  on  the  curvature  of  the 
pressure  distribution,  but  on  the  distribution  of  the  acceleration.  As  observed  in  figure  5,  the  pressure  and 
velocity  are  not  always  in  phase,  and  therefore  the  inclusion  of  the  this  term  can  induce  an  oscillation  in 
those  regions  where  the  fluid  acceleration  is  large,  such  as  the  arteries  and  veins. 

At  1.68s,  the  high-frequency  oscillation  of  the  arterial  pressure  begins,  which  was  also  observed  near  the 
beginning  of  the  heartbeat.  This  is  due  to  the  high  sensitivity  of  both  the  mitral  and  aortic  valves,  which 
open  or  close  at  each  timestep  because  the  pressures  in  the  heart  chambers  and  the  aorta  are  roughly  equal. 
If  the  valves  are  designed  to  only  close  when  the  downstream  pressure  is  more  than  7.5  mmHg  larger  than 
that  upsteam  of  the  valve,  this  chattering  of  the  arterial  pressure  signal  is  eliminated,  as  shown  in  figure  7(a). 
The  pressure  profiles  for  the  two  different  valve  sensitivities  are  qualitatively  very  similar.  Using  the  less 
sensitive  valves  that  close  at  the  higher  adverse  pressure,  there  is  a  period  of  increased  venous  pressure  from 
1.2s  to  1.34s.  This  is  a  result  of  the  venous  semi-lunar  valve  downstream  of  the  sensor  location  not  activating 
while  the  atria  contract.  In  fact,  using  the  less  sensitive  valve  model,  the  venous  semi- lunar  valves  are  not 
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(a)  (b) 

Figure  7:  Time  resolved  plots  of  (a)  pressure,  area,  and  (b)  velocity  at  three  locations  in  the  vasculature. 
Results  are  from  a  simulation  using  the  acceleration  model  (G  =  —du/dt).  In  this  model,  valves  only  close 
when  the  downstream  pressure  is  more  that  7.5  mmHg  larger  than  the  upstream  pressure. 


activated  during  the  entire  heartbeat. 

Comparing  the  results  of  the  two  different  acceleration  models,  a  signihcant  difference  exists  in  the 
magnitudes  of  the  arterial  pressure  profiles.  The  peak  arterial  pressure  using  the  zero-acceleration  model 
was  54  mmHg,  while  including  the  acceleration  model  more  than  doubles  it  to  147  mmHg.  This  is  attributed 
to  the  relatively  high  fluid  acceleration  in  the  arteries,  which  amplifies  the  maximum  pressure  at  1.33s  and  the 
minimum  pressure  at  1.4s.  Alternatively,  in  the  capillaries  where  there  is  relatively  little  fluid  acceleration, 
the  peak  pressure  only  increases  from  17.8  mmHg  in  the  zero-acceleration  case  to  22.1  mmHg  here.  In  the 
veins,  the  acceleration  model  decreases  the  minimum  pressure  to  35  mmHg  below  atmospheric,  whereas  it 
was  only  1.5  mmHg  below  atmospheric  when  using  the  zero-acceleration  model. 

The  velocity  profiles  obtained  using  the  G  =  —du/dt  acceleration  model  differ  from  the  zero  acceleration 
case  similarly  to  the  way  the  pressure  profiles  do,  as  seen  in  figure  6(b).  At  the  beginning  of  the  heartbeat, 
the  velocities  at  all  three  data  sensor  locations  are  non-zero  and  still  oscillating  from  the  previous  heartbeat. 
The  arterial  velocity  begins  to  increase  more  significantly  at  1.6s,  shortly  after  the  atria  begin  to  contract. 
At  1.25s,  when  the  ventricles  begin  to  contract,  the  arterial  velocity  increases  at  an  even  greater  rate  and 
the  high  frequency  oscillation  attributed  to  the  mitral  and  aortic  valve  sensitivities  is  eliminated.  At  this 
point  in  time,  the  velocity  in  the  capillaries  also  begins  to  increase  slowly.  The  arterial  velocity  peaks  at 
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1.35s,  and  the  capillary  velocity  peaks  later,  at  1.41s.  When  the  ventricles  begin  to  expand  at  1.45s,  the 
oscillating  pressure  observed  in  figure  6(a)  creates  a  period  of  low-amplitude  reversed  flow  in  the  arteries 
and  the  capillaries. 

The  venous  velocity  begins  to  increase  at  1.4s  when  the  atria  begin  to  expand.  It  peaks  at  1.6s,  and 
shortly  after,  at  1.63s,  the  capillary  velocity  peaks  a  second  time.  This  double  peak  in  capillary  velocity  is 
similar  to  that  the  zero-acceleration  case.  Due  to  the  oscillating  pressure  that  occurs  in  the  arteries  near 
the  end  of  the  heartbeat,  there  is  an  additional  period  of  positive  arterial  velocity  from  1.69s  to  2s,  during 
which  the  valve  sensitivity  creates  the  high  frequency  oscillation  in  the  arterial  velocity  profile.  As  seen  in 
hgure  7(b),  this  chattering  is  again  eliminated  by  using  valves  of  lower  sensitivity,  which  are  not  activated 
until  the  downstream  pressure  exceeds  the  upstream  pressure  by  more  than  7.5  mmHg.  Because  the  semi¬ 
lunar  valves  are  never  activated  in  this  case,  as  previously  mentioned,  a  period  of  reversed  flow  occurs  in  the 
veins  from  1.2s  until  1.37. 

In  table  2,  mean  and  maximum  velocities  at  the  three  sensor  data  locations  are  also  included  for  the 
acceleration  model  case.  As  expected,  inclusion  of  the  acceleration  model  increases  the  mean  and  maximum 
velocities  at  all  three  locations,  just  as  it  was  observed  to  increase  the  oscillating  pressure  magnitudes.  The 
mean  venous  velocity  is  slightly  greater  than  the  arterial  velocity,  similar  to  the  zero-acceleration  case,  but 
the  ratio  of  arterial  to  venous  maximum  velocities,  however,  is  lower  (261/173  =  1.5). 

5.3  Effects  of  system  elasticity 

After  establishing  the  baseline  distributions  of  pressure,  area,  and  velocity  through  the  model  vasculature, 
we  can  begin  to  explore  the  effects  of  varying  input  parameters  and  applying  additional  forcing.  In  hgure  8, 
we  show  the  time-resolved  pressure  in  the  arteries  and  capillaries  for  three  values  of  rigidity,  simulated  using 
both  the  zero  acceleration  and  the  (G  =  —du/dt)  model.  An  increase  in  rigidity  is  applied  as  a  decrease 
in  the  elasticity  (e)  in  the  equation  of  state  throughout  the  vasculature.  Using  the  G  =  0  model,  a  25% 
decrease  in  elasticity  cause  a  6%  increase  in  peak  pressure  in  the  arteries  and  a  15%  increase  in  peak  capillary 
pressure.  Halving  the  elasticity  caused  an  19%  increase  in  peak  arterial  pressure  and  a  58%  increase  in  the 
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Figure  8:  Effects  of  vasculature  rigidity  (a)  G  =  0  and  (b)  G  =  —du/dt 


G  =  0 

100%  e  75%  e  50%  e 

G  =  —du/dt 

100%  e  75%  e  50%  e 

Arteries 

54.6  57.8  (+6%)  62.9  (+15%) 

147.3  160.3  (+9%)  174.2  (+18%) 

Capillaries 

17.8  21.2  (+19%)  28.1  (+58%) 

22.1  25.7  (+16%)  32.1  (+45%) 

Table  3:  Peak  pressures  (in  mmHg)  for  different  elasticities  in  the  two  models,  and  the  percentage  increase 
over  the  baseline  value. 


capillaries.  Similar  increases  were  also  seen  when  using  the  G  0  model,  and  these  are  summarized  in 
table  3.  A  25%  decrease  in  elasticity  cause  a  9%  increase  in  peak  pressure  in  the  arteries  and  a  16%  increase 
in  peak  capillary  pressure,  while  a  50%  decrease  in  elasticity  caused  an  18%  increase  in  peak  arterial  pressure 
and  a  45%  increase  in  the  capillaries. 

When  including  the  acceleration  model,  another  effect  of  the  increased  rigidity  is  shown  in  in  figure  8(b). 
As  rigidity  increases,  not  only  do  the  peak  pressures  increase  over  the  baseline  case,  they  also  occur  slightly 
earlier.  The  loss  of  flexibility  advances  the  dynamic  response  of  the  system.  Furthermore,  from  1.37s 
to  1.46s,  the  increased  rigidity  appears  to  mitigate  some  of  the  heart  valve  chattering  effects  because  the 
arterial  pressure  profile  is  smoother.  Future  work  will  investigate  the  effects  of  varying  the  elasticity  along 
the  vasculature  and  the  system  sensitivity. 
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5.4  Effects  of  distributed  muscle  contractions 


Distributed  muscle  contractions,  as  described  in  §3,  were  applied  to  simulations  using  both  acceleration 
models,  and  the  results  were  evaluated  by  analyzing  time-averaged  flow  velocities  in  the  veins.  It  is  at  this 
location  where  the  physical  benefit  of  these  contractions,  represented  as  an  velocity  increase,  are  considered 
important  to  help  drive  blood  back  towards  the  heart  despite  the  relatively  low  venous  pressures. 

Test  cases  were  run  for  muscle  contractions  at  frequencies  in  the  range  0.5  Hz  <  /^  <  2.0  Hz  and 
maximum  amplitudes  in  the  range  0.005  <  ^AmjA  <  0.02.  A  summary  of  the  muscle  contraction  effects 
are  plotted  in  figure  9,  where  mean  venous  blood  flow  velocities  are  shown  for  four  contraction  amplitudes 
at  four  contraction  frequencies.  These  mean  velocities  were  averaged  over  six  heartbearts  to  ensure  than  an 
integer  number  of  muscle  contractions  were  included  in  the  timeseries  for  all  frequencies.  In  both  acceleration 
model  cases,  increasing  amplitude  increases  the  blood  flow  benefit  at  all  but  one  contraction  frequency.  At  I 
Hz,  there  is  no  benefit  to  venous  flow  according  to  the  zero-acceleration  model,  and  when  using  the  nonzero- 
acceleration  model,  muscle  contractions  at  this  frequency  actually  inhibit  venous  flow.  This  is  assumed  to 
be  a  result  of  contracting  in  phase  with  the  heart  pumping,  effectively  working  against  it. 

For  the  cases  studied  here,  the  lowest  frequency  of  muscle  contractions  (0.5  Hz)  created  the  greatest 
venous  flow  benefit  for  both  acceleration  models  at  almost  all  amplitudes.  Future  work  will  include  a  phase 
shift  between  the  muscle  contractions  and  heart  pumping  in  order  to  better  understand  how  to  optimize 
venous  blood  flow. 


6  Discussion 

A  new  low-dimensional  numerical  model  of  the  human  circulatory  system  is  introduced,  consisting  of  a  one¬ 
dimensional  axisymmetric  circular  elastic  pipe  with  periodic  boundary  conditions.  Computational  runtimes 
were  on  the  order  of  100  times  faster  than  realtime.  Five  liters  is  considered  a  normal  volume  for  an  adult 
human  circulatory  system  (Guyton  &  Hall,  2000),  and  the  model  presented  here  has  a  pipe  length  of  5m  and 
a  cross-sectional  area  distribution  that  yields  a  total  circuit  volume  of  4.7  L.  The  fluid  forcing  is  accomplished 
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Figure  9:  Mean  venous  flow  velocities  for  a  range  of  contraction  amplitudes  and  frequencies,  (a)  C?  =  0  and 
(b)  G  =  —dujdt 


through  an  effective  contraction  and  relaxation  of  the  local  area  near  the  heart  and  is  aided  by  a  set  of  model 
valves  in  the  heart  and  along  the  vascular  tree  to  prevent  expected  backflow. 

This  circulatory  system  model  performs  reasonably  well  and  predicts  some  large  scale  dynamic  quantities. 
The  magnitude  of  pressure  oscillations  is  approximately  correct  and,  as  expected,  it  is  largest  in  the  arteries 
and  decreases  precipitously  in  the  capillary  region.  Large  pressure  oscillations  were  also  observed  in  the 
veins  due  to  the  low  pressures  created  by  the  left  atrium  expansion.  Mean  velocities  at  points  in  the  arteries, 
capillaries,  and  the  veins  compared  well  with  reference  data  from  Charm  &  Kurland  (1974)  (table  2).  One 
main  discrepancy  is  the  over-prediction  of  mean  velocity  in  the  veins,  which  results  from  the  large  low- 
pressure  oscillations.  Future  model  development  that  resolves  the  correct  venous  pressure  signature  should 
also  decrease  the  local  mean  velocity. 

The  model  was  also  applied  using  a  modified  set  of  equations  that  included  the  pressure  drop  caused  by 
a  fluid  acceleration  term.  When  the  acceleration  effects  are  included,  the  time-resolved  pressures  exhibit  a 
larger  amplitude  in  the  the  arteries  and  veins,  where  the  fluid  acceleration  is  expected  to  be  high.  Also,  in 
these  regions  the  pressure  oscillates  at  a  frequency  that  is  higher  than  the  heartbeat  due  to  the  nonlinear 
feedback  of  the  fluid  acceleration  into  the  model.  This  formulation  also  makes  it  more  difficult  to  predict 
the  effects  that  improvements  to  the  model  geometry  and  characteristics  of  the  system  model  will  have  on 
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the  nonzero-acceleration  simulations. 


When  using  the  nonzero-acceleration  model,  a  very  high-frequency  pressure  oscillation  also  occurs  near 
the  beginning  and  end  of  the  heartbeat,  when  the  fluid  is  decelerating  and  the  pressure  in  and  around 
the  heart  chambers  is  fairly  constant.  The  chattering  of  the  pressure  (and  velocity)  comes  from  the  rapid 
opening  and  closing  of  the  mitral  and  aortic  valves.  If  the  valves  are  calibrated  to  only  close  at  a  larger 
adverse  pressure  gradient,  this  chattering  is  eliminated.  Future  work  will  refine  the  implementation  of  these 
valves  so  that  they  perform  more  smoothly  and  realistically. 

Changes  to  the  model  parameters  and  inputs  from  the  baseline  simulation  produce  realistic  responses  in 
the  circulatory  model.  An  increase  in  vascular  rigidity  by  25%  increases  the  arterial  pressure  maximum  in 
the  G  =  0  (G  7^  0)  simulation  by  6%  (9%).  Doubling  the  rigidity  increased  the  maximum  arterial  pressure 
by  15%  (18%).  The  application  of  distributed  muscle  contractions,  used  to  model  the  added  pumping  from 
skeletal  muscle  contractions,  was  shown  to  increase  the  mean  venous  velocities  in  the  model.  At  the  largest 
distributed  muscle  contraction,  with  maximum  amplitude  of  2%,  there  is  an  increase  of  20%  in  the  mean 
venous  velocity  using  the  zero-acceleration  model,  and  5%  using  the  nonzero-acceleration  model. 

There  are  several  geometric  and  mechanical  characteristics  of  the  current  model  that  need  to  be  refined 
as  development  continues,  despite  the  current  qualitative  similarities  on  the  macro  length  and  time  scales. 
The  relative  cross-sectional  areas  and  lengths  of  the  arteries,  capillaries,  and  veins  are  currently  defined 
heuristically,  so  that  the  total  volume  and  length  are  physiologically  realistic.  Also,  the  two  halves  of  the 
circulatory  loop  (systemic  and  pulmonary)  are  symmetric,  even  though  it  is  known  that  this  is  not  the  case 
in  human  physiology.  The  physical  pulmonary  circuit  is  considerably  shorter  and  offers  less  resistance  to  the 
output  of  the  right  ventricle  than  the  systemic  circuit  offers  to  the  left  ventricle.  Future  development  of  the 
model  will  focus  on  the  realistic  calibration  of  the  quantities  such  as  segment  length,  area,  distribution  of 
the  weighted  viscosity  and  elasticity,  and  valve  representation. 

With  respect  to  the  model  of  the  actuated  heart  chambers,  it  is  not  in  the  scope  of  this  work  to 
correctly  capture  the  complex  three-dimensional  dynamics.  Rather,  we  will  use  calibrated  phenomenologies 
to  capture  the  results  of  experiments  and  more  detailed  simulations.  For  this  model,  it  is  important  to 
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recreate  a  physiologically  realistic  cardiac  output,  and  future  work  to  refine  the  model  of  heart  pumping 
will  focus  on  replicating  the  quantity  of  fluid  pumped,  but  not  necessarily  on  geometrically  matching  the 
physical  heart. 

A  primary  concern  about  the  large-scale  dynamic  behavior  of  the  model  is  that  it  does  not  yet  correctly 
capture  a  high  mean  pressure  in  the  arteries.  In  reality,  the  high  pressure  created  by  the  contraction  of  the 
heart  keeps  the  aorta  and  the  large  arteries  at  a  pressure  near  100  mmHg  above  atmospheric.  We  have  so 
far  attempted  to  address  this  issue  in  two  ways.  First,  the  local  equilibrium  area  was  set  to  be  greater  in  the 
large  arteries,  which  initiated  a  larger  arterial  pressure,  but  this  relaxes  to  the  lower  equilibrium  pressure 
within  two  heartbeats.  Consequently,  the  local  pressure  is  much  lower  than  the  equilibrium  pressure  in  the 
arteries,  causing  a  significant  decrease  in  the  local  area. 

A  second  approach  to  correct  this  discrepancy  was  to  increase  the  weighted  viscosity  in  the  capillary 
region.  This  not  only  drives  up  the  pressure  in  the  arteries  as  the  heart  pumps,  but  at  high  enough  values, 
the  arterial  pressure  relaxes  to  equilibrium  at  a  rate  slower  than  the  heart  rate.  The  arterial  pressure  then 
oscillates  around  a  higher  mean,  as  desired.  This  approach,  however,  significantly  decreases  the  low  pressure 
on  the  venous  end  of  the  vasculature.  The  enhanced  resistance  in  the  capillaries  achieved  by  a  higher  viscosity 
does  increase  the  arterial  pressure,  but  at  the  cost  of  unrealistically  low  venous  pressure  and  high  venous 
velocities.  Further  research,  augmented  with  additional  expertise,  is  needed  to  resolve  this  issue. 

In  addition  to  the  continued  improvement  of  this  circulation  model,  other  future  work  is  planned  to 
expand  the  scope  of  this  project.  Currently,  a  similar  formulation  is  being  used  to  develop  a  model  of  the 
respiratory  system.  These  two  independent  models  will  then  be  coupled  together  through  a  model  of  the 
respiration  process,  involving  the  transport  of  gases  (O2  and  CO2)  across  the  blood  gas  barrier.  We  also  plan 
to  develop  an  independent  fluid  dynamic  model  of  the  lymphatic  system,  although  a  new  formulation  may 
be  needed  to  appropriately  capture  its  diffusion-dominated  dynamics.  The  lymph  (interstitial  fluid)  comes 
from  a  slow  “leak”  of  material  out  of  circulatory  system,  and  the  lymphatic  system  returns  it  at  the  venous 
ends  of  the  systemic  vasculature.  This  process  provides  another  natural  coupling  between  the  circulatory 
and  lymphatic  system  models. 
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The  refinement  of  the  circulatory  system  model,  the  development  of  independent  respiratory  and  lym¬ 
phatic  models,  and  the  careful  coupling  among  them  will  occur  simultaneously.  It  is  as  yet  too  soon  to  state 
that  the  interaction  of  the  relatively  simple  models  will  be  able  to  recreate,  on  some  level,  the  intricate  and 
nonlinear  behavior  of  the  human  body  as  a  whole,  but  the  progress  is  promising. 
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